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Abstract 

The use of quartic force fields (QFFs) has been shown to be one of the most effective 
ways to efficiently compute vibrational frequencies for small molecules. In this paper we 
outline and discuss how the simple-internal or bond-length bond-angle (BLBA) coordi- 
nates can be transformed into Morse-cosine(-sine) coordinates which produce potential 
energy surfaces from QFFs that possess proper limiting behavior and can effectively de- 
scribe the vibrational (or rovibrational) energy levels of an arbitrary molecular system. 
We investigate parameter scaling in the Morse coordinate, symmetry considerations, and 
examples of transformed QFFs making use of the MULTIMODE, TROVE, and VTET 
variational vibrational methods. Cases are referenced where variational computations 
coupled with transformed QFFs produce accuracies compared to experiment for funda- 
mental frequencies on the order of 5 cm -1 and often as good as 1 cm -1 . 

Keywords: quartic force fields, morse-cosine coordinates, variational vibrational 
methods, vibrational configuration interaction theory, vibrational frequencies 


1. Introduction 

The quartic force field (QFF) is a simple and elegant means by which one can ade- 
quately describe the near equilibrium potential surface for a molecular system. A QFF is 
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simply the fourth-order Taylor series approximation of the anharmonic potential of the 
form: 

^ = 2 FijAiAj + - FijkAiAjAk + — FijkiAiAj AfcA; (1) 

ij ijk ijkl 

with displacements Aj and force constants F. tJ . Methods utilizing QFFs have been 
shown to provide accuracies often as good as 1 cm -1 with experiment and routinely as 
close as 5 cm -1 in describing the fundamental vibrational frequencies for various systems 
[1-18] though QFFs based on coupled cluster singles, doubles, and perturbative triples 
[CCSD(T)] [19] energies produce such accuracies partly due to a cancellation of errors in 
the electronic structure method, as has been discussed [20-22] . Vibrational perturbation 
theory (VPT) [23-25] at second-order (VPT2) is most often combined with QFFs for 
determination of rovibrational properties. Much success is reported in the literature 
for reproducing experimental frequencies and reporting previously unobserved modes 
[26-33]. These studies and many others showcase how beneficial molecular frequencies 
computed with QFFs and VPT2 are to the chemistry and physics communities. 

There have been some cases, however, that highlight where VPT2 can fail, and these 
differ from known weaknesses in QFFs. When VPT2 is not a valid choice of vibrational 
method, other tools must be employed to solve the nuclear Schrodinger equation. For 
example, Martin and Taylor [34] found that the strongly anharmonic N— H stretching 
mode in the H 2 NN molecule cannot be adequately described with VPT2 but can be using 
variational procedures. This was consistent with an earlier study that found a similar 
N— H stretch anharmonicity for the HNO molecule [3]. Further work [35] analyzed the 
role that variational computations can play in improving the description of anharmonic 
frequencies and found that non-negligible anharmonic effects in hydrogen peroxide and, 
to a lesser degree, hydrogen persulfide are not properly accounted for within VPT2. 
Hence, expansion-style wavefunction-based vibrational methods are sometimes the only 
means of reliably describing vibrational frequencies beyond the harmonic approximation. 

Variational methods for the computation of anharmonic vibrational frequencies can 
take many forms, and many are connectivity-specific or, in the very least, formulated 
for a specific number of atoms. Work on tetra-atomic systems increased the amount 
of chemistry that can be accurately modeled by quantum chemical procedures. In the 
VTET approach by Schwenke [36] and the DVR(n) methods by Mladenovic [37], explicit 
terms for the kinetic and potential energy terms are derived based on how the atoms 
interact with one another within a chosen or developed set of vibrational basis functions. 
Recent work [38-41] has shown how these approaches are beneficial in resolving issues 
with computed vibrational frequencies. 

Other methods and codes, especially those of the newest generation, are designed to 
be more general and are built around terms that can be applied to systems of an arbi- 
trary number of atoms without regard for the specific connectivity or symmetry. The 
TROVE program [42] constructs a general vibrational Schrodinger equation and expands 
the kinetic and potential energy operator terms in the Hamiltonian with polynomials. 
An active space of vibrational states is chosen to contain the expansion, and the potential 
energy function (which can be input from QFFs) for certain geometrically defined co- 
ordinates completes the Hamiltonian. Diagonalization of the Hamiltonian results in the 
vibration-rotation energies and wavefunctions, where the latter are the descriptions of 
the actual molecular behavior. This scheme also gives purely rotational states and mixed 
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vibration-rotation states. TROVE has been utilized since its development to provide line 
lists for common systems like NH 3 [43, 44] and even to resolve some of the anharmonicity 
issues experienced with peroxide-like systems such as HSOH [45, 46]. 

Another, generally-applicable program designed to treat anharmonic vibrational com- 
putations is MULTIMODE (MM) [47-50] . This program makes explicit use of vibrational 
configuration interaction (VCI) theory in order to describe the rovibrational wavefunc- 
tion. Within the Hamiltonian, the potential is defined as a hierarchical combination 
of intrinsic potential terms, V^ N \Qi) for coordinates Qi. N is the number of modes 
considered in a single potential term, and the maximum of N leads to the so-called 
mode-representation (MR) of the total VCI potential. The coupled potential terms pro- 
duce the VCI potential: 


J p is the total angular momentum of a given cardinal direction ( x , y , or z) denoted by p 
or a in Eq. 3; 7r p is the total vibrational angular momentum of the same direction; fi pa 
is the inverse of the moment of inertia tensor for the given geometric coordinates; and 
Q is the set of all normal coordinates. With the potential and Hamiltonian defined, the 
vibrational self-consistent field (VSCF) equations can be formulated and solved. These, 
in turn, are fed into the VCI equations. Both sets of equations and their solution steps are 
discussed in Ref. [48] . In VCI, full-CI is achieved when N or, equivalently, AMR equals 
the total number of modes the molecule may possess. For a four-atom system, full-CI is 
6MR from 3 N — 6. Table 1 gives a comparison for MM and TROVE for ammonia where, 
at the full computational level for each method, the two agree to better than 0.4 cm -1 
for every frequency. 

However, all of these post-harmonic, non-VPT2, variational methods require the in- 
put of a potential function. QFFs often reliably describe the potential near equilibrium. 
Global or semi-global surfaces are useful, as well, but they are much more computa- 
tionally costly in the formulation of the potential. However, the larger expansions of the 
variational wavefunctions routinely require states where bond length-bond angle (BLBA) 
or simple-internal coordinates are known to break-down when QFFs are used to define 
the potential. Hence, in order for these methods to be effective, they must have what is 
called “correct limiting behavior.” 

2. Transforming to a Morse-Cosine Coordinate System 

2.1. Correct Limiting Behavior 

In simple-internal or BLBA coordinates, QFFs are known to produce non-physical 
results when attempting to describe behavior at distances just beyond the equilibrium 
geometry and at values larger than the QFF displacements [52-54]. This is especially 



( 2 ) 


The above VCI potential is then included in the Watson Hamiltonian [51, 24]: 



( 3 ) 
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true for the bond lengths where the energy, depending on the case, can either grow con- 
tinually toward positive infinity as the bond length is stretched [3], or turn over and 
deteriorate below the minimum and continue toward negative infinity as other formula- 
tions of coordinates have been shown to do [55-57]. For the case of the O— H stretch in 
water, shown in Fig. 1, the BLBA coordinate QFF, denoted by the black curve, exhibits 
the unchecked increase toward positive infinity. This is also present in the N— H and C— H 
stretches for ammonia and methane in Figs. 2 and 3. Hence, the resulting variational 
wavefunction becomes incapable of describing behavior for distances anything more than 
slightly larger than the equilibrium bond lengths. If the reference geometry is not close 
to the actual potential minimum or if the artificial barrier for a turnover to negative 
infinity is not tall or wide enough, this problem can be more pronounced, depending on 
the variational method. 

For bond stretches, the Morse potential converges to an asymptotic limit where the 
atoms become noninteracting. In order for a potential energy surface (PES) to possess 
“correct limiting behavior,” the energies of the displacements for longer bond lengths 
must converge to some value in a similar way. This does not have to be the actual Morse 
potential limit even though such would definitely be the desired case, but the PES must 
have a similar shape. The incorrect limiting behavior of the PES given by a QFF along a 
bond length coordinate can be easily rectified by transforming the QFF into a coordinate 
system in which correct limiting behavior is conserved. 

An easily understood transformation of coordinates to address PES and QFF correct 
limiting behavior was brought forth by Simons, Parr, and Finlan (SPF) [57] originally 
for diatomic molecules. In these coordinates, R, the equilibrium bond length, r e , is used 
to modify the displaced bond length, r, in the form: 



These SPF coordinates often prohibit the potential from turning over towards negative 
infinity as the related Dunham [55, 56] or simple BLBA coordinates are known to do. 
In SPF coordinates, r replaces r e as the denominator from their Dunham counterparts. 
SPF coordinates do not directly model the Morse-like behavior at longer bond lengths 
[57]. In fact, they often diverge noticeably from the Morse-potential shortly beyond the 
equilibrium geometry and are not guaranteed to converge to any limit. The difference 
between SPF coordinate-based energies and the Morse limits is great enough that at 
large enough bond lengths, highly non-physical behaviors are modeled, but this is not as 
great as is often observed in BLBA coordinates that may approach positive infinity much 
faster [57, 3]. Hence, SPF coordinates often will not give correct limiting behavior but do 
so in quite a different manner as compared to Dunham or BLBA representations. That 
being said, SPF coordinates can force the weighting of the wavefunction towards the 
actual bottom of the well overcoming one issue with coordinates without correct limiting 
behavior. As a result, computations that require longer bond lengths in the potential 
are possible with SPF coordinates, but they are not as robust as would be desired. 

2.2. Form of the Morse Coordinate 

Other formulations of coordinates attempt to rectify the issues that SPF cannot. 
Watson [51] proposed a means of describing bond lengths with Morse-like behavior even 
before SPF coordinates were introduced. However, it was not until Meyer, Botschwina, 
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and Burton [58] popularized this approach, followed shortly by Carter and Handy [59], 
that it began to gain acceptance. These new coordinates, R, are Morse-like and are 
defined [58] as: 


R = 


1 — e 




ft 


( 5 ) 


where ft is a parameter that can be optimized for the system of interest. It is this style 
of bond length coordinate that has been applied most straightfowardly to QFFs for use 
in describing the anharmonic potential for variational computations. 

Later, a different but related formulation of the Morse-like coordinates was suggested 
by Dateo, Lee, and Schwenke [3] : 


R = 1 — e 


r (r— r e ) 


( 6 ) 


where a r is again an optimizable parameter that can vary for each bond length present 
in the computation. The general parameter, a;, was defined as: 


OLi 



( 7 ) 


where these diagonal elements of the three- and two-body force constant matrices are 
still in internal coordinates and have not yet been transformed. Such a definition of the 
a parameter forces the transformed force constant, F ,; lt , to be zero. Unlike Eq. 5, Eq. 6 
does not include the parameter as the denominator term in its formulation, but this still 
leaves the equation with convergence to some limit and is nearly the exact form as the 
Morse potential itself, mostly by design. Even so, coordinates defined either by Eq. 5 or 
Eq. 6 will behave similarly though they converge to slightly different limits. Hence, both 
have correct limiting behavior at least in a semi-quantitative sense. The red curves in 
Figs. 1-3 showcase how Morse-type coordinates (as defined by Eq. 6) differ from BLBA 
coordinates, especially as the bond length is increased. The inflection points as the red 
curves trend toward a limit can be observed. As a result, Morse coordinates for bond 
stretching modes are essential when using QFFs together with variational anharmonic 
vibrational calculations. 


2.3. Cosine for Simple Bending Angle Coordinates 

It is not just the bond lengths that need modification, however. Methods like VCI in 
MULTIMODE allow the modes to couple together within the potential and are further 
susceptible to inaccuracies in the PES brought about by the descriptions of the modes 
from the chosen coordinate type. Potential curves along a single bond angle must be 
periodic in nature. However, BLBA coordinates may overshoot the vibrational energy at 
various points on the potential, especially for values close to n [3]. These non-physical 
descriptions using BLBA coordinates can also be overcome with transformations which 
need not be exceptionally complicated. The standard bond angle displacement coordi- 
nate, termed the theta coordinate, is defined below: 

Q = 6-e e . ( 8 ) 

Coupled with Morse coordinates for the stretches, this coordinate system is denoted 
Morse-theta. 
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However, as indicated, it is often necessary to have the potential exhibit periodic 
behavior, and this may be achieved by taking the cosine of each angle: 

0 = cos(0) — cos (6 e ). (9) 

The use of Morse coordinates and these cosine transformed bond angle coordinates 
as a set is referred to as Morse-cosine coordinates. QFFs in this coordinate system 
have been previously shown to produce highly accurate frequencies when coupled with 
variational computations [3, 6, 9, 10, 12-18]. As an example, the QFF-based fundamental 
frequencies for ammonia have been computed with the VTET program using both Morse- 
cosine and Morse-theta coordinates. These are given in Table 1. For the bond stretching 
frequencies, the change in the bend coordinate has an effect of a negligible 0.15 cm -1 or 
less, but the change between Morse-cosine and Morse-theta coordinates is a noticeable 
1.7 cm^ 1 or more for the bending modes. Furthermore, Fig. 4 shows the difference in 
behavior between the theta bond angle coordinate in black and the cosine coordinate in 
red for water. For this triatomic molecule, the major difference in the two is at n where 
the theta coordinate overshoots the energy and can introduce errors in the vibrational 
calculation. 

Lastly, close examination of the results in Table 1 shows that even though only a 
QFF was used for ammonia, the variational treatment for solving the nuclear Schrodinger 
equation leads to a symmetric double-well potential, and hence the inversion splitting 
of all energy levels, which is a hallmark for the ammonia molecule. (The reader should 
note that for the E modes, we report different parities due to the inversion splitting, and 
hence they are not identical.) In fact, comparison of the various splitting values between 
the QFF results, whether from TROVE or VTET, with the values obtained from a 
highly accurate global PES, HSL2 [39, 40], shows that they are in remarkably good 
agreement. The inversion splitting values obtained using the MULTIMODE program 
are still reasonably good, but not as accurate as obtained with VTET and TROVE. In 
this case, the ammonia QFF used in the variational calculations has an inversion barrier 
of 1778 cm -1 versus 1785 cm -1 for the HSL2 global PES. 

2.f. The Choice of Alpha 

Since the a value in the Morse coordinate of Eq. 6 is an arbitrary parameter, this 
value can be calibrated to give different results. As such, a can be chosen to exactly 
describe the dissociation limit determined from high-accuracy theoretical or experimental 
data. For cases where such data is not known or cannot be determined, a standard and 
accurate definition of a is valuable. Dateo, Lee, and Schwenke [3] proposed the definition 
of a given above in Eq. 7, but this choice of a has not been systematically studied for 
its performance in computations of the vibrational frequencies with Morse coordinates. 
Tables 2-4 highlight how changes in the a parameter affect the frequency computed for 
the fundamentals of NH 3 , H 2 0, and CH 4 , respectively. 

One obvious choice of a would be the one that gives the proper dissociation energy 
( D e ). The computed D e for each QFF of a given molecule, as listed in Tables 2-4, is the 
asymptotic limit that the potential surface approaches as r increases toward positive in- 
finity. These D e values are computed from relaxed optimizations which means they allow 
the other bonds to respond and reach a local minimum energy geometry as the principle 
bond is stretched. When scaling the computed D e with regards to the experimental D e , 
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it is shown that the standard (Eq. 7) choice of a performs quite well in describing the 
dissociation energy. A ±5% range for the Eq. 7 choice in a contains the experimental 
D e for our three test cases. The results for a are within 3.5% of the experimental values 
for ammonia, 1.84% for water, and 2.45% for methane. Water would need a smaller a 
to match the experimental D e while ammonia and methane would need a larger a value 
to do so. 

In Table 2, the ammonia fundamental vibrational frequencies are computed with 
TROVE. The QFF used is fitted from energies determined from within the HSL2 global 
PES [39, 40] at individual displacements of 0.005 A for the bond lengths and 0.005 
radians for the bond angles and is transformed into the Morse-theta coordinate system. 
The fundamental vibrational frequencies computed are within 4.5 cm -1 or better with 
the spectroscopically accurate HSL2 results. Furthermore, a variance of ±5% from the 
standard a value computed by way of Eq. 7 changes the stretching frequencies by only 
1.24 cm -1 and 1.61 cm -1 . The degenerate deformation bends are affected by 0.07 cm -1 
or less. These results indicate that the value of a is well-behaved for the definition given 
in Eq. 7. Compared with the HSL2 results for the fundamentals in Ref. [39, 40], many of 
the modes, such as V 2 and 1/3 , could benefit from a decrease in the a parameter. Others, 
like ^4 ± compare nicely for the unsealed a in the “1.0” column of Table 2. Regardless, 
this choice of a created without any assumption or experimental data input is able to 
give correct limiting behavior whose asymptotic limit is very close to the experimental 
D ei and also is an integral part in the generation of fundamental frequencies to fairly 
high accuracy. 

The behavior of the frequencies for water and methane with MM VCI using Morse- 
cosine coordinates in Tables 3 and 4 is similar to the ammonia results computed with 
TROVE. However, the fundamentals for water are computed with the cC QFF from 
Ref. [9]. The latter QFF is so named here with the convention utilized later [13, 18] 
since it is defined from CCSD(T) [19] (described above) utilizing aug-cc-pCVTZ, aug- 
cc-pCVQZ, and aug-cc-pCV5Z basis sets [60-62] extrapolated to the complete basis set 
(CBS) limit with a three-point formula [63]. The so-called LMT QFF used for the 
methane computations is from Ref. [4] and is defined from CCSD(T)/cc-pVQZ energies 
at points defined from displacements like those in ammonia. 

From Tables 3 and 4 it can be seen that the fundamentals for water differ from 
experiment by less than 2.5 cm -1 and differ by less than 4 cm -1 for methane. The 
variance between frequencies for values of a within the ±5% range for water and methane 
are on the order of 0.1 cm -1 for the angle bends and less than 2.0 cm -1 for the bond 
stretches. While ammonia and water have many modes that would benefit slightly from 
a decrease in the choice of a, many modes of methane would perform better with an 
increase in a. Further, there are several modes for these three test cases where the 
frequencies compare favorably with experiment for the Eq. 7 choice of a. Hence, the 
current choice of a, which does not require any prior knowledge of the system analyzed, 
can give reliable fundamental vibrational frequencies provided the QFF is also accurate. 

The differences in the a values are visually depicted in Figs. 1-3 where the blue, 
red, and green curves correspond to 0.7, 1.0 (unsealed), and 1.3 a scalings for the X— H 
stretches in water, ammonia, and methane, respectively. The three different colored 
curves in each figure begin to diverge a few tenths of an Angstrom beyond the equilibrium 
bond length. The blue 0.7 a potential curves are higher in energy than the standard red 
curve beyond the divergence point, and the 1.3 a green curves are lower in energy. This 
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behavior is expected from the negative exponential definition of the Morse coordinate 
in Eq. 6. Not only does the choice of a defined in Eq. 7 readily produce accurate 
fundamental frequencies, it also yields results for the D e values that are very close to the 
experimental quantity and, as a result, correct limiting behavior for systems for which 
experimental data is not available. 

2.5. More Complicated Structures 

There are occasions where simple Morse-cosine transformed coordinates for a QFF 
are not adequate to treat a system. Two examples are planar, cyclic structures and linear 
structures. In cases where the equilibrium bond angle coordinate value is close to 0 or 
7r, sine coordinates of the form: 


0 = sin(0) — sin(0 e ) (10) 

are necessary. In the fitting of the force constants, derivatives of the coordinates appear 
in denominator terms. At the equilibrium geometry of a linear structure, for instance, 
^ cos(7r) = — sin(7r) = 0. Hence, a non-zero gradient term is necessary to avoid an 
infinity value arising from a 0 in the denominator. Usage of these coordinates in con- 
junction with cosine coordinates for the other angles with bonds that are greater than 0 
and less than 7r leads to the so-called Morse-cosine-sine coordinates. The usage of these 
coordinates has been shown to be valuable in linear Good molecules like CCH" [10] where 
the sine coordinate is necessary to accurately treat the n bending mode. 

Furthermore, the high levels of redundancy for some coordinates, like the C— C bonds 
within benzene together with the C— C— C bends, can create havoc when manipulating 
QFFs. As a result, symmetry adaptation of the simple internal coordinates must be 
used in highly symmetric structures before construction of the QFF and transformation 
to Morse-cosine(-sine) coordinates. Symmetry adaptation has been utilized for highly 
symmetric structures for some time now [1, 4, 64-66]. Further, recent work on c-C3H3 + 
[12] exemplifies the need for symmetry-adaptation in the coordinate system for a Morse- 
type coordinate, as well. The initial QFF is constructed in terms of symmetry-adapted 
internal coordinates making full use of the D$h point group symmetry. For example, 
complete treatment of the C— C bond lengths becomes a single term, S i = (Ri + i? 2 + 
Rs)/V 3, since there is a redundancy between the three C— C bond lengths and the three 
C— C— C bends. A full listing of these coordinates is given in Ref. [12], and all are linear 
combinations of the redundant set of bond lengths, bond angles, and out-of-plane bends, 
where the latter are defined as the motion of one hydrogen atom relative to the carbon 
ring. Usually, the symmetry internal coordinate QFF is transformed directly into the 
Morse-cosine coordinate system using the definitions in Eqs. 6 and 9. However, in the 
case of C-C3H3+, it is necessary to use a symmetry adapted set of Morse-cosine-sine 
coordinates due to the redundancy between the C— C bonds and the C— C— C bends. In 
this case, one can derive formulae to directly transform the symmetry internal coordinate 
QFF into a symmetry adapted Morse-cosine-sine coordinate QFF, or what was done in 
Ref. [12] was to re-fit the QFF using the original displacement energies represented in the 
symmetry adapted Morse-cosine-sine coordinate system. In this way, the new symmetry 
adapted Morse-cosine-sine QFF has proper limiting behavior, and it also has exact D^h 
symmetry. 


3. Degenerate Vibrational Modes 


Generally, VPT2 computations can treat degenerate vibrational modes exactly as long 
as the QFF is properly constructed in terms of the symmetry coordinates, and this has 
been done many times ( e.g . Refs. [1, 4]). However, for variational vibrational calculations 
the proper treatment of degenerate vibrational frequencies is more complicated. Different 
variational vibrational methods treat symmetry in different ways. This is mainly a 
response to the formulation of the method itself and how the symmetry relationships fit 
within the equations. For those methods that are general, especially with regards to the 
number of possible atoms included in the computation, the symmetry treatment must 
be generalized in some fashion as similar to electronic structure computations [67]. 

However, neither variational vibrational nor electronic structure codes are regularly 
coded to treat symmetries of point groups beyond D 2 /j and its subgroups. D 2 h and its 
subgroups all possess purely one-dimensional irreducible representations (irreps) which 
give straightforward mathematical relationships to be coded. Hence, degenerate vibra- 
tional modes may not be properly treated. As a result, linear combinations of lower 
symmetry irreps must be employed to describe the degenerate modes of higher symme- 
tries. In the MULTIMODE VCI program, C 2 V and its subgroups, C s and C 2 , can be 
treated. Higher-order point group molecules can also be included, but they must use the 
largest Abelian subgroup available. As a result of this lack of explicit symmetry treat- 
ment, there are three different ways in which inadequate treatment of degenerate modes 
in higher-order symmetries can be introduced: the Hamiltonian, the wavefunction, and 
the basis size, or, in other words, inadequate convergence thresholds. 

3.1. The Hamiltonian 

In order to maintain symmetry, the Hamiltonian must be exact or approximated in 
such a way that components of a degenerate vibrational mode are treated in the same 
way [24, 68]. An example of where this problem was identified in an electronic structure 
method is given in Ref. [69] for the OPT2, open-shell perturbation theory, method. For 
the vibrational problem, when considering systems of high symmetry, there can exist 
errors in the formulation of the kinetic energy terms for point groups with multidimen- 
sional irreps if care is not taken when approximating the kinetic energy operator. Most 
notably, however, the QFF (the potential term in the Hamiltonian) may be easily forced 
to properly treat the symmetry relationships of the system. This was shown in a few early 
cases [1, 4, 64], and QFF symmetry relationships have now been published for molecules 
of the form X 3 {D 3h ) [65], X 4 (T d ) [65], XV 3 (<7*0 [70, 1, 27], XY 3 (D 3h ) [71, 7], XY 4 
(T d ) [4, 64], XYZ (C 0 ov ) [72, 73, 2], XY 2 (D„ h ) [72, 71], X 2 V 2 (fl^) [74, 72, 75], and 
x 3 y 3 (D 3h ) [12], 


3.2. The Wavefunction 

The second consideration of symmetry is from the wavefunction perspective. The 
symmetry modes for MULTIMODE are defined by the user in the input. The MM VCI 
potential (Eq. 2) must be totally symmetric since the Hamiltonian, by definition, is totally 
symmetric [76], and the QFF terms must be also. Hence, the different irrep-labeled mode 
interrelation to mimic higher-symmetry is already defined when the potential (Eq. 2) is 
created. This forces the wavefunctions themselves to describe the desired symmetry 
states. The TROVE program operates in a similar fashion. The explicit degeneracies 
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present in the higher symmetry can be lost when the irreps are transformed to lower 
symmetry. Degenerate states are comprised of linear combinations of various irreps. As 
a result, certain pieces of the necessary linear combination may be discarded if those 
states are above a certain inclusion threshold. Hence, more terms in the VCI expansions 
and more basis functions per irrep are typically necessary to fully describe the actual 
degeneracies. Within a 3, 4, or 5MR VCI computation, a degenerate pair of vibrational 
modes may not be fully contained within the wavefunction since all combinations of three, 
four, and five modes are included in the highest substitution level of the wavefunction 
terms where other terms, including the degenerate partner, may be excluded. Hence, 
higher substitution levels in the Cl wavefunction allow for a more complete description 
of the vibrational wavefunction. At the full Cl level, all of the possible excited states in 
the wavefunction are included and degenerate levels will be described properly, but only 
a few interesting molecules can be computed at this level. 

Table 5 highlights how larger numbers of terms in the wavefunction improve the vari- 
ational computation of the frequencies for degenerate modes. For the 4MR v 3 frequencies 
of ammonia computed with MULTIMODE, the and values differ by 23.92 cm -1 . 
However, simply increasing the number of virtual states included in the wavefunction by 
using 5MR decreases this difference by two orders of magnitude to 0.59 cm -1 . 6MR is 
of course full Cl for ammonia, and the differences between the 5MR and 6MR results 
are small. As discussed previously, the E modes for ammonia do not come out exactly 
degenerate because different parities of the inversion split energy levels are reported. 

The degenerate modes of methane, given in Table 6, show a better case for the use 
of larger expansions in the VCI wavefunction. The doubly degenerate 02 frequencies 
differ by 0.20 cm -1 for 3MR, 0.02 cm -1 for 4MR, and are degenerate to better than 
0.01 cm -1 for 5MR. The z/4 degenerate frequencies behave similarly as the components 
for V 2 - The 1/3 frequencies differ by only 0.13 cm -1 for the 5MR computations, while 
the 3MR frequencies are in disagreement by 8.95 cm -1 . BH3 (given in Table 7) also 
shows improvement for the higher substitution level in the VCI wavefunction. The QFF 
defined from the CBS extrapolated energy refined for non-Born-Oppenheimer and scalar 
relativisitic effects (CBR) has already been employed with VPT2 [7] with high accuracy. 
Using the same QFF, the degenerate E' mode components of BH 3 are nearly exactly 
degenerate at the 5MR level, while the other E' mode, ^4, varies between the computed 
energy levels by less than 0.01 cm -1 . Hence, increasing the size of the wavefunction will 
improve the description of the degenerate modes such that discrepancies between the 
two or more computed components in a system with a single-well minimum can agree to 
better than 0.01 cm -1 . 

3.3. Convergence Thresholds 

The issue with energy-level degeneracies computed with lower symmetries also re- 
quires the use of very tight convergence criteria in the vibrational wavefunction. This 
approach has been proven effective for electronic structure computations [69]. Tighter 
convergence thresholds must also be applied to the construction of the wavefunctions 
in the rovibrational computations. Work on the cyclic (D 3 h) form of C3H 3 + [12] has 
showcased the reliability in using QFFs and C2„-based VCI to describe the rovibrational 
properties of a highly symmetric system. BH3 in Table 7 also shows how larger basis 
sets can improve the agreement in degenerate components, where the 04 components of 
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1196.6889 and 1196.6936 cm -1 shrink the 0.0053 cm -1 gap to 0.0012 cm -1 with fre- 
quencies of 1196.6059 cm -1 and 1196.6071 cm -1 . This demonstrates that the size of 
the underlying basis set is important for convergence of computed energy levels and 
maintaining the symmetry of different components of a degenerate mode. 

4. Combination Bands, Overtones, and Higher-Energy States 

As previously mentioned, variational vibrational approaches making use of Morse- 
cosine QFFs have been able to produce accuracies for fundamental vibrational frequencies 
of routinely 5 cm -1 and as close as 1 cm -1 in some cases [3, 6, 9, 10, 12-18]. However, 
comparative data for the overtones and combination bands is not as common. The 
overtones and combination bands are often much less intense and of less significance in 
the analysis of experimental vibrational spectra. However, they can play a defining role 
in the spectra of certain systems, especially smaller molecules with fewer fundamental 
vibrational frequencies. 

If correct limiting behavior and proper coordinates are necessary for the fundamental 
vibrational frequencies, they are more important for the combination bands and over- 
tones. This should be apparent since the energy levels for combination and overtone 
bands will be higher up the potential well. However, examination of the VPT2 term 
value expression shows that even harmonic frequency errors are magnified for higher en- 
ergy levels [24, 25], and this will be apparent in variational calculations of vibrational 
energy levels, as well. In the following section we compare the transition energy for VCI 
combination and overtone bands relative to experiment for the water molecule (using the 
aforementioned Morse-cosine cC QFF), followed by a comparison of the VCI and VPT2 
results. 

4-1. Comparison of VCI with Experiment 

As an example, several theoretical (cC 3MR VCI MULTIMODE) and experimental 
(Ref. [77]) fundamental vibrational frequencies, overtones, and combination bands of 
water are listed in Table 8. Since water only has three degrees-of-freedom, 3MR is full- 
CI. The cC VCI fundamentals differ from the corresponding experimental frequencies by 
an average of 4.08 cm -1 . As a result, some of this difference gets compounded in the 
overtones, especially for the v 2 bending frequency. The state with the most difference 
in the computed result as compared to experiment is the 4649.42 cm -1 3^2 overtone 
which is 17.37 cm -1 lower than the 4666.790 cm -1 experimental result. The next-largest 
experimental difference is 8.31 cm -1 from 2iz 2 , which is an overtone of the same mode. On 
the opposite end of the error spectrum, the tq+tq combination band exhibits the smallest 
theoretical difference from experiment, 1.85 cm -1 . These examples appear to indicate 
some level of compounding or cancellation of the error as will occur with combination 
and overtone bands 

The VCI overtones and combination bands will compound or cancel the errors some, 
but only to a certain degree as appears in the above cases. However, whether the error 
compounds or is canceled to some extent for a higher-lying vibrational state is difficult to 
predict in advance since errors in the harmonic part of the calculation may be different 
(and the opposite sign) to those for the anharmonic part. To illustrate, even though the 
computed 3661.86 cm -1 iq differs from experiment by 4.81 cm -1 , the computed 7205.20 
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cm -1 2iq differs from experiment by only 3.66 cm -1 , and 3iq differs by an even less 3.30 
cm -1 . Moreover, iq + has an error of 2.15 cnr -1 from experiment. This is less than 
half of what the error is for either of the stretching fundamentals. In both of these cases, 
errors for the combination and overtone bands were smaller than for the fundamentals, 
showing a cancelation of errors for the higher-lying states. 

The average difference from experiment for the overtones included in Table 8 is only 
4.97 cm -1 . Furthermore, besides those frequencies involving overtones of u 2 , no state 
differs from experiment by more than 5.0 cm -1 with an average error of 3.12 cm -1 . Thus, 
the VCI results in Table 8 show that a QFF transformed into a Morse-cosine coordinate 
system can yield accurate vibrational frequencies for combination and overtone bands, 
though the results for 2v 2 and, especially, 3v 2 suggest that caution should be exercised 
in evaluating the results. 

4-2. Comparison of VCI and Experiment with 2 nd Order Perturbation Theory 

For tightly bound molecules, most fundamental vibrational frequencies show strong 
correlation between VPT2 computed with SPECTRO [78] and VCI/QFF computed with 
MULTIMODE with differences in the frequencies most often less than 5 cm -1 . Even so, 
there are many exceptions to this rule as has been shown in the work on the conformers of 
the HOCO radical [13, 14] where the VPT2 and VCI torsional modes differ by more than 
25 cm -1 . (It has been shown that VPT2 is fortuitously more reliable in its computation of 
this mode as compared to DVR(6) [41].) However, for tightly-bound molecules, examples 
like this almost always are a result of the molecule possessing a large-amplitude motion, 
and that is the case here. Similar situations have been reported for NH 3 , HNO, and 
H 2 NN, as discussed earlier. 

In looking at Table 8, it would appear that VPT2 performs slightly better in the 
description of the fundamental vibrational frequencies than VCI where both methods 
make use of the cC QFF [9]. VPT2 is nearly exact for is 2 , but the 1591.96 cm -1 VCI 
frequency is only 2.82 cm -1 less than experiment. VPT2 reports the ZPE to be 4662.72 
cm -1 while VCI is much lower at 4644.68 cm -1 , nearly a factor of four closer to the 
experimentally inferred 4638.39 cm -1 ZPE. For iq, VCI and VPT2 are within 1.0 cm -1 
of each other with VPT2 only 0.67 cm -1 higher. As a result, VCI is slightly closer 
to the experimental iq frequency at 3657.04 cm . Finally, VPT2 and VCI are nearly 
coincident for their computation of the i/ 3 antisymmetric stretch, and are less than 5.0 
cm -1 higher than experiment. The average error compared to experiment for the three 
fundamentals of water computed with VPT2 is 3.56 cm -1 , less than 4.08 cm -1 with VCI. 
Closer inspection reveals that VPT2 does a slightly better job describing the bending 
frequency than VCI, but does not perform as well as VCI for the two stretching modes. 

This trend continues for the overtones, combination bands, and higher-energy states. 
The v 2 mode and its overtones are better described with VPT2 than they are with VCI, 
but the difference between VPT2 and experiment is consistently greater than the dif- 
ference between VCI and experiment for all of the other modes. The average error for 
VPT2 for all of the states given in Table 8, including the fundamentals, is 12.14 cm -1 , 
7.27 cm -1 higher than that from VCI. The maximum VPT2 difference from experiment 
is 38.32 cm -1 above the experimental 8761.582 cm -1 2u\ + v 2 frequency while the next- 
largest error is 36.29 cm -1 for 2iq. Hence, it is essential to use methods beyond VPT2 
when computing overtones, combination bands, and higher-energy states, and VCI cou- 
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pled with a QFF transformed into a Morse-cosine coordinate system appears to be a 
viable alternative. 

5. Conclusions 

Procedures to utilize QFFs in MULTIMODE or other variational approaches like 
those found in TROVE or VTET have been under-appreciated even though they have 
existed for many years [51] and have been utilized for decades [58, 59]. In this work 
we have shown how QFFs can be used effectively in the computation of vibrational 
frequencies even with variational methods. 

In order to achieve correct limiting behavior in QFFs, the coordinates must be trans- 
formed into some scheme that properly describes the energy beyond the vicinity of the 
potential energy minimum. For bond lengths, this is a modified form of the Morse poten- 
tial which allows for asymptotic convergence to some dissociation limit. The limit varies 
with the value of the scaling parameter, a, but the experimental result is contained with 
a ±5% deviation of the a value defined in a general form nearly twenty years ago [3] 
(Eq. 6). Additionally, a ±5% deviation from the standard alpha is shown to contain the 
exact experimental frequency for a given fundamental from our sample set of ammonia, 
water, and methane. The bond angles also need to be transformed in order to describe 
properly the potential at values approaching and surpassing 7 r. This transformation is 
especially important when the bond angle, or a similar coordinate like a torsion, displays 
some periodic behavior. The most effective means to treat bond angles is through the 
use of a cosine coordinate (Ecp 9), but some cases, such as high symmetries or specific 
out-of-plane bends, require the use of the related sine coordinate (Eq. 10). 

Symmetry considerations are often necessary in defining the QFF and utilizing them 
with variational methods. For high symmetry cases {i.e., where there are degenerate 
vibrational frequencies), the potential (i.e., the QFF), even in the Morse-cosine(-sine) 
coordinate system, can be forced to obey symmetry constraints such that appropriate 
degeneracies are exact. Such symmetry relationships amongst the QFF force constants 
have been derived and published for a whole host of molecules. Thus, the main source 
of loss of exact degeneracy in a variational vibrational calculation is due either to the 
variational method treating the components of a degenerate mode differently (as can 
happen with VCI if less than a full Cl is used) or simply due to a lack of convergence. 
In both cases, these sources of error can be overcome with a large enough calculation. 

When using the same QFF, variational methods like VCI are often coincident with 
VPT2 in computing the fundamental vibrational frequencies. A few exceptions have 
been noted where VPT2 performs better than VCI [41], but this has only been shown as 
a better correlation between VPT2 and DVR(6) , which is explicitly coded to treat tetra- 
atomics. VCI often produces frequencies just as close or closer to experiment than VPT2 
as shown here and previously [12, 13, 16, 17]. For overtones and combination bands, 
VCI generally provides a more complete description of the state wavefunction and and a 
more accurate frequency than the simple perturbational approach in VPT2. For water, 
the average difference from experiment in the VCI computation of the fundamentals, 
overtones, and combination bands is less than 5 cm' 1 while the average error for VPT2 
is more than 12 cm -1 . 

In summary, the use of QFFs in conjunction with variational computations is a valid 
means by which one may study vibrational or rovibrational frequencies provided the QFF 
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is transformed into an appropriate coordinate system. QFFs require relatively few points 
to create a potential surface while sacrificing little in terms of accuracy. This reduces 
the total computational cost while still producing meaningful results. With the growth 
of generalized variational vibrational codes, the use of Morse-cosine transformed QFFs 
allows for various systems to be examined that would be too costly otherwise. 
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Figure 1: The 1-dimensional potential energy curves of the unrelaxed O— H stretch in water: Black is 
BLBA, blue is Morse-cosine with the 0.7 scale factor for a, red is Morse-cosine without a scaling, and 
green is Morse-cosine for 1.3 a. The black horizontal line is the O— H symmetric stretch fundamental for 
BLBA coordinates while red, slightly below, is the same but with unsealed a Morse-cosine coordinates. 



Table 1: Comparison between vibrational methods and coordinate systems for the fundamentals (in 
cm -1 ) of NH 3 . a 


mode 

Morse-theta 

VTET 

Morse-cosine 

VTET 

Morse-cosine 
MM 6MR 

Morse-cosine 

TROVE 

HSL2 6 

VTET 

Non- adiabatic 
HSL2 VTET 


4 

3332.98 

3333.09 

3332.07 

3332.32 

3336.48 

3336.10 

K 


3334.06 

3334.13 

3333.19 

3333.38 

3337.47 

3337.08 

A'i 

4 

927.10 

928.80 

928.94 

928.65 

932.64 

932.44 

K 

V 2 

964.31 

965.67 

965.32 

965.45 

968.43 

968.16 

K 

0“ 

0.829 

0.82 

1.14 

0.82 

0.795 

0.793 

E' 

4 

3440.86 

3441.00 

3440.21 

3440.18 

3444.04 

3443.63 

E" 

v -i 

3441.21 

3441.35 

3440.92 

3440.55 

3444.40 

3443.98 

E' 

4 

1623.75 

1625.67 

1625.20 

1625.33 

1626.75 

1626.28 

E" 

V A 

1624.89 

1626.80 

1626.70 

1626.46 

1627.85 

1627.37 


“Unless otherwise noted, all computations employ the QFF components of the HSL2 
global PES from Refs. 39, 40. 

6 These computations were executed with the HSL2 global PES from Ref. 39, 40. 
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Figure 2: The 1-dimensional potential energy curves of the unrelaxed N— H stretch for ammonia: Black 
is, again, BLBA, blue is Morse-theta with the 0.7 ck, red is Morse-cosine without a scaling, and green is 
Morse-cosine for 1.3 a. 
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Figure 3: The 1-dimensional potential energy curves of a single unrelaxed C— H stretch in methane: 
Black is BLBA, blue is 0.7 a Morse-cosine, red is non-scaled Morse-cosine, and green is Morse-cosine for 
1.3 a. 
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Figure 4: The 1-dimensional potential energy curves of the unrelaxed H— O— H bond angle for water: 
Black is BLBA, and red is from Morse-cosine coordinates. 
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Table 5: Comparison of the NH 3 4MR, 5MR, and 6 MR MM fundamental vibrational frequencies (in 
cm -1 ). 


Symmetry 

Mode 

4MR 

5MR 

6MR 

A 


3331.95 

3332.06 

3332.07 

K 

v i 

3332.99 

3333.18 

3333.19 

A[ 

A 

928.13 

929.06 

928.94 

K 


963.98 

965.47 

965.32 

K 

0“ 

1.13 

1.14 

1.14 

E' 

V t 

3439.20 

3440.35 

3440.21 

E" 

^3 

3463.12 

3440.94 

3440.93 

E' 

-4 + 

1619.74 

1624.98 

1625.20 

E" 

f 4 

1620.91 

1626.48 

1626.70 


Table 6 : Comparison of the CH 4 3MR, 4MR, and 5MR MM fundamental vibrational frequencies (in 
cm -1 ). 


Symmetry 

Mode 

3MR 

4MR 

5MR 

Tli 

v 1 

2915.08 

2916.28 

2916.41 

E 


1532.56 

1532.50 

1532.63 



1532.76 

1532.52 

1532.63 

t 2 


3018.87 

3018.29 

3017.97 



3027.82 

3017.49 

3018.10 



3027.82 

3017.49 

3018.10 

t 2 

f 4 

1312.69 

1312.91 

1313.04 



1312.66 

1312.94 

1313.04 



1312.66 

1312.94 

1313.04 


Table 7: Comparison of the CBR a QFF BH 3 3MR, 4MR, 5MR, and 5MR-l b MM fundamental vibrational 
frequencies along with those from VPT2 C and experiment (in cm -1 ). 


Symmetry 

Mode 

3MR 

4MR 

5 MR 

5MR-P 

VPT2“ 

Expt. d 

A 1 

v 1 

2500.4673 

2500.9143 

2500.9318 

2500.8695 

2502.3 


Tli 

F2 

1145.5123 

1145.5366 

1145.5444 

1145.5445 

1147.2 

1147.4986 

E' 


2602.8146 

2602.8933 

2602.9349 

2602.8806 

2602.1 

2601.5743 



2606.8095 

2602.8471 

2602.9348 

2602.8771 



E' 

f 4 

1196.6687 

1196.6686 

1196.6889 

1196.6059 

1196.5 

1196.66 



1196.5841 

1196.6734 

1196.6936 

1196.6071 




“The CBR QFF is from Ref. 7 and is the CBS extrapolated, Born-Oppenheimer cor- 
rected, scalar relativistic QFF. 

^The “-1” indicates that these 5MR results are using a basis set with more than twice as 
many functions. 

“The VPT2 computations are from Ref. 7. 

d The experimental v 2 and f 4 frequencies are from Ref. 83 with i/ 3 from Ref. 84. 
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Table 8: H 2 O fundamentals and first few overtones and combination bands from cC VPT2, cC 3MR 
(full-CI) VCI, experiment, and the differences from experiment for VPT2 and VCI (all in cm -1 ). 


Mode 

VPT2 

VCI 

Exp.“ 

Exp. - VPT2 

Exp. - VCI 

ZPE 

4662.72 

4644.68 

4638.39 

24.33 

6.29 

^2 

1594.54 

1591.96 

1594.746 

0.21 

2.79 

2^ 

3153.87 

3143.32 

3151.630 

-2.24 

8.31 

Vl 

3662.53 

3661.86 

3657.053 

-5.48 

-4.81 

"3 

3760.92 

3760.58 

3755.929 

-4.99 

-4.65 

3v 2 

4677.98 

4649.42 

4666.790 

-11.19 

17.37 

v 1 + Vi 

5240.84 

5236.83 

5234.978 

-5.86 

-1.85 

v 2 + V3 

5335.68 

5334.04 

5331.265 

-4.41 

-2.77 

Vl + 2^2 

6783.92 

6770.35 

6775.093 

-8.83 

4.74 

2^2 + ^3 

6875.22 

6868.94 

6871.520 

-3.70 

2.58 

2v\ 

7237.83 

7205.20 

7201.540 

-36.29 

-3.66 

v 1 + V3 

7254.98 

7251.97 

7249.818 

-5.16 

-2.15 

2^3 

7423.49 

7453.23 

7445.045 

21.56 

-8.18 

2v\ + v 2 

8799.90 

8763.69 

8761.582 

-38.32 

-2.11 

V\ + v 2 + V3 

- 

8809.05 

8806.999 

- 

-2.05 

v 2 + 2v 3 

8978.46 

9008.41 

9000.136 

21.68 

-8.27 

3^i 

- 

10596.39 

10599.686 

- 

3.30 


a Reference 77. 
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